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ABSTRACT 

The numerical investigation of Bondi-Hoyle accretion onto a moving black hole has a long 
history, both in Newtonian and in general-relativistic physics. By performing new two- 
dimensional and general-relativistic simulations onto a rotating black hole, we point out a 
novel feature, namely, that quasi-periodic oscillations (QPOs) are naturally produced in the 
shock cone that develops in the downstream part of the flow. Because the shock cone in the 
downstream part of the flow acts as a cavity trapping pressure perturbations, modes with fre- 
quencies in the integer ratios 2 : 1 and 3 : 1 are easily produced. The frequencies of these 
modes depend on the black-hole spin and on the properties of the flow, and scale linearly with 
the inverse of the black-hole mass. Our results may be relevant for explaining the detection of 
QPOs in Sagittarius A*, once such detection is confirmed by further observations. Finally, we 
report on the development of the flip-flop instability, which can affect the shock cone under 
suitable conditions; such an instability has been discussed before in Newtonian simulations 
but was never found in a relativistic regime. 
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1 INTRODUCTION 

Non-spherical accretion flows are astrophysically relevant in all 
those situations where strong winds with a small amount of angu- 
lar momentum are able to transport considerable amounts of matter 
towards an accreting compact object, e.g., a black hole or a neu- 
tron star. Examples include the case of massive X-ray binaries in 
which a compact object accretes from the wind of an early-type 
star, or the case of young stellar systems orbiting in the gravita- 
tional potential of their birth cluster and accreting from the dense 
molecular interstellar medium. One of the most celebrated and 
best st udied types of non-spher i cal accretion flows is the Bondi- 
Hoyle dHovle & Lvttletonll 19391 ; iBondi & Hovldl 19441) which, we 
recall, develops when a black hole moves relative to a uniform 
gas cloud. The Bondi-Hoyle flow has bee n the subject of sev- 
eral numerical in vestigations, start i ng fro m lMatsuda et al.l fl987l) 



region of the flow and it manifests in the oscillation of the shock 
cone from one side to the other of the accretor. 

In the relativistic regime, the first two-dimensi onal simula 
tions of Bondi-Hoyle accreti on were performed by iPetrich et al 



1989) and subsequently by iFont & Ibanezl dl998bT ): iFont et al 



1993. 1 19991) . These investigations considered flows both in ax- 



Frvxell & Taarnl Jl988h ; ISawada et alTdl989h 



and followed by I 

iBenensohn et all ( fl997n T As a sum mary of the bulk of w ork done 
so far on this topic, Table 1 of lFoglizzo"et ID d2005h reports 
an overview of published numerical simulations of Bondi-Hoyle- 
Lyttleton accretion over the last 30 years and lists more than 40 
works. While providing detailed descriptions of the morphology of 
the gas capture in the supersonic regime, several of the above men- 
tioned investigations were aimed at verifying the occurrence of the 
so called flip-flop instability. This consists of an instability to tan- 
gential velocities of the shock cone that forms in the downstream 
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isymmetry and in the equatorial plane; interestingly none of these 
works showed the occurrence of the flip-flop instability, so that no 
evidence existed, prior to this work, about the development of the 
instability also in a curved spacetime. More recently, a fully gen- 
eral relativisti c investigation of th e Bondi-Hoyle accretion has been 
considered bv lFarris et al .1 ( |2010| ) in the context of the merger of su- 
permassive black-hole binaries; in this case, however, no discussion 
of the flip-flop instability was presented. 

However, the occurrence of the instability is not the only rel- 
evant physical process that may manifest in Bondi-Hoyle accre- 
tion flows. An additional one, that does not seem to have been 
considered so far, is related to the possibility that the shock cone 
traps oscillation modes, thus producing Quasi Periodic Oscillations 
(QPOs). The intuition is indeed simple yet rather suggestive. The 
typical shocked cone of the Bondi-Hoyle accretion is likely to pro- 
vide a natural cavity for the development and confinement of os- 
cillation modes of sonic nature. If the accretor is a compact object, 
Newtonian hydrodynamics is a good approximation only at large 
radial distances, while it can lead to misleading conclusions when 
studying the flow evolution close to the event horizon. In this paper 
we perform two dimensional general relativistic hydrodynamics 
simulations by focusing on the dynamical behavior of the Bondi- 
Hoyle shock cone under a wide range of parameters. The QPOs 
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that we have found have typical frequencies in the range from 
10" 5 Hz to 10" 3 Hz for a black hole with M 

l2 1 



10 6 M Q (and in 



the range from 1 Hz to 10 Hz for a black hole with M = 10 M ). 
As a result, they can become relevant for the interpretation of QPOs 
observed both in the galactic center and in high-mass X-ray bina- 
ries. Our analysis has also shown that the flip-flop instability does 
occur even in the relativistic framework, and we have investigated 
how such instability can interplay with the QPO phenomenon, by 
suppressing or exciting specific modes of oscillations. 

The plan of the paper is the following: in Section [2] we pro- 
vide an overview of the hydrodynamical equations, of the numeri- 
cal methods adopted and of the physical set up of the problem. In 
Sec. [3] we describe our results, while Sec. [4] contains a discussion 
of the potential applicability to two specific astrophysical cases. Fi- 
nally, Sec. [5] is devoted to the conclusions of our work. We use a 
geometrized system of units with G = c = 1. 



2 MATHEMATICAL FORMULATION 

2.1 General relativistic hydrodynamics equations 

We study Bondi-Hoyle accretion of a perfect fluid in the curved 
background spacetime of a rotating black hole. The energy mo- 
mentum tensor of the fluid has the usual form 



Table 1. Initial models adopted in numerical simulation. From left to right 
the columns report: the name of the model, the black hole spin a, the asymp- 
totic flow velocity Voo, the asymptotic Mach number .Moo and the inner 
boundary r; n of the grid. The final simulation time is set equal to 10000M 
in all of the models. The reference adiabatic index adopted is T = 4/3, 
though we have also considered different values to explore the effect on the 
results. 



Model 


a/M 


Voo 


Moo 


r in (M) 


Al 


0.9 


0.001 


0.01 


1.78 


A2 


0.9 


0.1 


1 


1.78 


A3 


0.9 


0.2 


2 


1.78 


A4 


0.9 


0.3 


3 


1.78 


A5 


0.9 


0.4 


1 


1.78 


A6 


0.9 


0.5 


5 


1.78 


A7 


0.9 


0.6 


6 


1.78 


B 


0.5 


0.4 


4 


2.1 


CI 


0.0 


0.001 


0.01 


2.1 


C2 


0.0 


0.1 


1 


2.1 


C3 


0.0 


0.2 


2 


2.1 


C4 


0.0 


0.3 


3 


2.1 


C5 


0.0 


0.4 


1 


2.1 


C6 


0.0 


0.5 


5 


2.1 


C7 


0.0 


0.6 


6 


2.1 



= hpu^u" + pg M 



(1) 



where u M is the four velocity of the fluid, h, p and p are the specific 
enthalpy, the rest-mass density and the pressure, respectively. All 
of these quantities are measured by an observer comoving with it M , 
while g^ v is the metric of the spacetime. By adopting the 3 + 1 
formalism of general relativity and after choosing Boyer-Lindquist 
coordinates (t, r, 0, (f>), the line element of the metric is written as 



ds 



where 

E 2 = 
A = 
A = 



2Mr\ 2 



— =5 — sm t)atd(f> 



r + a 2 cos (9) , 



? sin 



r - 2Mr - 
(r + a ) 



2 . .2 
a A sm i 



(2) 
(3) 
(4) 



with a and M being the spin and the mass of the black hole. The 
lapse function and the shift vector of the metric are given, re- 
spectively, by ol = (E 2 A/A) 1/2 and ft = (0, 0, -2Mar/A). 
When solving numerically the general-relativistic hydrodynamic 
equat ions it is importa nt to write them in a conservation 
form jBanvuls et alii 19971) 



<9U dF 1 



(5) 



where U, F l and S are the vectors of the conserved variables, of 
the fluxes and of the sources, respectively, while q 1 is the gen- 
eralized coordinate in the i— th direction. When performing two- 
dimensional numerical simulations in the equatorial plane, Eq.(|5) 
reduces to 



dU dF r OF 4 

dt dr dd> 



(6) 



The conservative variables, written in terms of the primitive vari- 
ables (p, v 1 ,p), are 



U 



D 












T 





L y/liphW 1 



^yWp 

^ P hW 2 Vr 

^phW 2 V^ 



Wp) 



(7) 



where V 1 = u l /W + /3 J /ct is a spatial vector whose indices 
are raised and lowered through the spatial metric jij and it rep- 
resents the three-velocity of the fluid with respect to the Eu- 
lerian observer associated to the 3 + 1 splitting of the metric. 
W = (1 - 7ijT / V j )" 1/2 is the Lorentz factor of the fluid and 
7 = det (7i,) = S 4 sin 2 8 /a 2 is the determinant of three metric. 
Although lFont et aT] d 19991) already reported fluxes and sources for 
the hydrodynamical equations in the Kerr metric, we explicitly list 
them here for convenience as 



F r = 



a(V r - /3 r /a)D 
a{(V r - p r /a)S r + V7P} 

a(V r - /3 r /a)S^ 
a{(V r ~p r /a)T + ^/V r p} 



(8) 



a{(V* - 



-pt/aQD 
-p*/a)S T 
- 0*/a)S* + 
P*/a)T + y/W*p} 



(9) 



and 



S = 







a^y(T rt d r a - a(T rt r^ + T"T* r + T^T%)) 

(10) 

The equation of state that we adopt is that of an ideal gas, 
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Figure 1. Rest-mass density distribution (left panel) and Lorentz factor distribution (right panel) for a black hole with spin a = 0.9 and a supersonic flow with 
asymptotic velocity Voo = 0.5, and at time t = 10000 M. In the downstream region a clear shock cone forms which is distorted by the non vanishing spin of 
the black hole. 



namely 



= (r-l)pe, 



(11) 



where F and e are the adiabatic index and the specific internal en- 
ergy, respectively. The reference adiabatic index used in the simu- 
lations is r = 4/3, though we have also considered different values 
to explore the effect on the results. 

We have solved the system of equations l[5} through high- 
resolution shock-capturing schemes (HRSC) based on approximate 
Riemann solvers. A minmod linear algorithm for the reconstruction 
of the left and right states at each interface between adjacent nu- 
merical cells is adopted, while the numerical fluxes are computed 
with the Marquina flux formula. More specific d etails about the nu- 
merical scheme can be found in lDonmea d2004l) . 

As a final remark we note that because we are not here consid- 
ering any contribution coming from the cooling of the shocked gas, 
the oscillations we will d iscuss are not rela t ed to th ose reported 
by iMolteni et al.l dl996h : IChakrabarti et al.l d2004h : bkuda et al.l 
d2007t) . In those studies, in fact, the presence of cooling processes 
is a necessary ingredient for the appearance of oscillations along 
the standing shocks. 



2.2 Initial Conditions, Boundary Conditions and 
Approximations 

We perform numerical simulations on the equatorial plane, i.e., 8 = 
7r/2. The initial velocity fie l d is g iven in terms of an asymptotic 
velocity Voo as in lFont et al.Ul99a) , i.e., 




7"- Voo COB((j>) , 

7^* Voo sin(0) . 



(12) 
(13) 



These relations guarantee that the velocity of the injected gas at 
infinity is parallel to the a;— direction, while V 2 = W l = 
everywhere in the flow. The value of Vx> can then be chosen to in- 



vestigate the different regimes of the flow and to consider therefore 
subsonic or supersonic accretion. 

During the evolution additional matter is injected supersoni- 
cally from the outer boundary in the upstream region with the same 
analytic prescription of dl2t and (13), thus reproducing a contin- 
uous wind at large distances. The initial density and pressure pro- 
files are adjusted to make the sound speed equal to a required value, 
which we choose to be c Sj0 o = 0.1. In practice, once c St00 is cho- 
sen, we set the density to be a constant (i.e., p = 1) and then 
the pressure is derived from the relativistic definition of the sound 
speed, i.e., p = c? )0O p(T - 1)/(T(T - 1) - c^T). A brief de- 
scription of the initial models is reported in Table [T] where models 
O refer to Schwarzschild black holes and different fluid injection 
velocities, while models A* and B refer to Kerr black holes respec- 
tively with spin a/M — 0.9 and a/M — 0.5, again spanning dif- 
ferent fluid injection velocities. 

The computational grid consists of N r x N$ uniformly spaced 
zones in the radial and angular directions, respectively, covering a 
computational domain extending from r m i n (reported in Table [TJ 
to r max = 43 M and from to 2-7T. For our fiducial simulation we 
have chosen N r = 512 and N$ = 256, but have also verified that 
the qualitative results (i.e., the appearance of the QPOs or of the 
instability) are not sensitive to the resolution used or to the location 
of the outer boundary, which has been moved to r max ~ 80 M in 
some tests. 

The boundary conditions in the radial direction are such that 
at the inner radial grid point we implement outflow boundary con- 
ditions by a simple zeroth-order extrapolation (i.e., a direct copy) 
of all variables. At the outer radial boundary, on the other hand, we 
must distinguish between the upstream region, with n /2 < 4> < 
3/2-7T, and the downstream region, with — n/2 < <f) < tt/2. In the 
upstream region we continuously inject matter with the initial ve- 
locity field, while in the downstream region we use outflow bound- 
ary conditions. Finally, periodic boundary conditions are adopted 
along the <f> direction. 
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Finally, we note that in order to validate the code, we have 
ca rried out a number of co mparisons with the ECHO code presented 
in Del Zanna et al ] ( l2007h . which, among other differences, uses 
a non-uniform grid structure. The results obtained from the two 
codes have not shown significant differences. 



3 NUMERICAL RESULTS 

3.1 Properties of the stationary pattern 

The main features of Bondi-Hoyle accretion were first investigated 
in the rela tivistic framework through a number of num e rical simu- 
lations by IPetrich et al.l Jl989h ; lFont & Ibafiej Jl998bl) : lFont et all 
dl 9981.1 19991) ." The overall results obtained can be summarized as 
follows. 

As already shown by several authors, when a homogeneous 
flow of matter moves non radially towards a compact object, a 
shock wave will form in the neighborhood of the accretor. Depend- 
ing then on the properties of the flow, namely on the adiabatic in- 
dex and on the asymptotic Mach number M x , the shock can either 
reach very close to the accretor or be at a certain distance from it 
[see, e.g., iFoglizzo et al .1 j2005h l. More specifically, for any given 
value of .Moo, there is a critical adiabatic index T ~ 2, below 
which a shock wave of conic shape, i.e., a "shock cone", forms in 
the downstream part of the flow, i.e., downstream of the accretor. 
Matter from the upstream region crosses continuously the shock 
front, it undergoes a strong deceleration and, if it is inside the ac- 
cretion radius, is ultimately accreted. As a result, the maximum 
rest-mass density in the downstream region is always larger than 
the corresponding one in the upstream region and, consequently, 
the mass accretion rate is significantly non-spherical and larger in 
the downstream part of the flow. Interestingly, when suitable coor- 
di nates are used (e.g ., Kerr-Schild coordinates), it has been shown 
bv lFont et al.l l ll999h that the shock even penetrates the event hori- 
zon. On the other hand, for increasingly stiff fluids and thus for val- 



ues of the adiabatic index larger than the critical one, a bow shock 
wave forms in the upstream part of the flow, i.e., upstream of the 
accretor. In this case the accretion is almost spherical and because 
the bow shock does not reach into the black-hole's horizon, it is 
often referred to as a "detached" shock. 

The shape of the shock cone depends sensitively on the spin 
of the black hole. If the black hole is not spinning, then the shock 
cone is perfectly symmetric about the = direction. On the other 
hand, if the black hole is spinning, the induced frame-dragging ef- 
fect produces a "wrapping" of the shock cone, as already shown 
and discussed in great detail by |Fontetal1 l ll999l) . 

Our simulations essentially confirm this picture, and the main 
features of the stationary pattern are summarized in Figure [T] 
which shows the two dimensional rest-mass density distribution 
(left panel) and the Lorentz factor distribution (right panel) for 
model A6 with a/M = 0.9 and Voo = 0.5. As evident from the 
iso-density contours, the accretion in the upstream region, is es- 
sentially spherical accretion, while a very well defined shock cone, 
characterized by a strong density gradient, forms in the downstream 
region. The shock cone, where the accretion rate is the largest, is 
clearly distorted by the spin of the black hole and its opening an- 
gle slightly decreases with increasing distance from the center. We 
have also found that spinning black holes have maximum rest-mass 
density in the shock cone which is higher than for non spinning 
black holes. The right panel of Fig. Q] on the other hand, shows that 
the overall flow is only very mildly relativistic, with a maximum 
Lorentz factor Wmmx < 3, even ahead of the shock front in the 
vicinity of the black hole. It is also interesting to note that accretion 
in the shock cone takes place with a very small Lorentz factor. Mat- 
ter inside the shock cone can accrete onto the black hole only for 
radii smaller than a given accretion radius r a , defined as the radius 
where the accretion rate (as a function of radius) becomes negative. 
The expression r a ~ 2M/(V^ + (? St0a ), which is valid in the case 
of a purely-spherical accretion, provides a good estimate of the ac- 
cretion radius in the Bondi-Hoyle accretion. As an example, for a 
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model with Voo = 0.5 and c Sj00 = 0.1, we found r a = 6.1, while 
the predicted one is r a = 7.6. 

To complement the information in Fig.[T] Fig. [2] shows one- 
dimensional profiles at time t = 10000 M of the rest-mass density 
at different radial shells for an asymptotic flow velocity Voo =0.5. 
More specifically, the left panel refers to Schwarzschild black hole, 
while the right one to a rotating black hole with spin a/M — 0.9. 
Clearly, a sharp transition in the density exists at the border of the 
shock cone and this is particularly evident for the nonrotating black 
hole. Note also the amount of distortion which is present in the case 
of a rotating black hole; although this is in part due to our choice 
of Boyer-Lindquist coordinates, similar distortions are present also 
when using different and bet ter suited coordin ates, such as Kerr- 
Schild [see the discussion in lFontetal.l i fl999l) l. As already men- 
tioned above, the angular location of the shock is only weakly de- 



pendent on the radial distance, with the size of the shock cone re- 
ducing only slightly when moving away from the black hole. As a 
final remark we note that the large degree of symmetry shown in the 
left panel of Fig.|2]provides a convincing evidence of the abilities 
of the code to maintain the symmetry across </> = 0. 

3.1.1 The flip-flop instability 

As shown through numerical simulations performed by several au- 
thors over the years in Newtonian physics, the Bondi-Hoyle accre- 
tion flow is subject to the so called flip-flop instability, namely an in- 
stability of the shock cone to tangential velocities, which manifests 
in the oscillation of the shock cone from one side of the accretor to 
the ot her. Such instability was first discovered by iMatsuda et al.1 
in 2-D Newtonian axisymmetric simulations, and later 
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confirmed and furt h er investigated b y Frvxell & Taaml 1 19881) 



were performed by 



Ishii et al 



ISawadaetal](ll989l). Benensohn et ; _. 

dl999h . I Pogorelov et al .1 J20001) . Thre e di mensiona l simu lations 



I b y IFryxell & Taaml H1988I). 
all i l 19971) . iFoglizzo & Ruffed 



i ll 9931) and iRuffertl Jl997h . who 
confirmed the occurrence of the instability, although with defor- 
mations of the shock cone only very close the accretor. In spite 
of all these investigations, however, the nature and the physical 
origin of the flip-flop in stability remain obscure. According to 
IFoglizzo & Rufferl dl999h . for example, local instabilities, such as 
the Rayleigh-Taylor or the Kelvin-Helmholtz instabilities, should 
not play a significant role and canno t acco unt for the flip-flop 
instability. On the other hand, Soker (1990) showed, through a 
WKB analysis, that the two-dimensional axisymmetric accretion 
flow (not the one considered in our work) can be unstable against 
tangential modes, as well as against radial modes. When extending 
these results to three-dimensional simulations, however, the persis- 
tence of such instabilities is not fully clarified. Overall, it is fair to 
say that, although the flip-flop instability has been confirmed nu- 
merically both in two and three-dimensional simulations, it is still 
unclear whether the physical mechanisms driving the instabilities 
in the two cases are the same or not. Fi nally, to fur t her co mplicate 
the scenario, no evidence was found bv lFont et al.l i 19981) of such 
instability within a relativistic framework, although the authors ar- 
gued that an instability might develop for very large values of the 
asymptotic Mach number. 

More recently, IFoglizzo et al. I <2005h performed a very de- 
tailed analysis of the unstable behavior of Bondi-Hoyle accretion 
flows, suggesting that the instability may be of advective-acoustic 
nature, both in the case of shocks attached or detached (bow 
shocks) to the accretor. Moreover, they also stressed that, though 
physical, the instability is likely to be triggered or dumped by nu- 
merical effects such as the carbuncle phenomenon at the shock, the 
boundary condition at the surface of the accretor and the grid size. 

Although a detailed analysis of the flip-flop instability is be- 
yond the scope of this paper, which is rather focused on the devel- 
opment of QPOs in Bondi-Hoyle accretion flows, our simulations 
confirm for the first time the occurrence of the instability even in 
a relativistic regime. Fig. [3] for instance, reports four snapshots of 
the rest-mass density in a long term evolution of a model when the 
instability is fully developed. The shock cone in the downstream 
region oscillates back and forth in the orbital plane in this represen- 
tative model with black hole spin a/M = 0.5 and Voo = 0.4. The 
reason why such unstable behaviour was not noticed bv lFont et al.l 
i 19991) may be due to the short term evolutio n that they were force d 
to consider at the time, or, as suggested by IFoglizzo et at] J2005b . 
because of the high values of 2M/r a = (VjJ, + ef i00 ) that they 
considered. 

Fig. [4] shows additional key features of the flip-flop instabil- 
ity. The top panel, in particular, reports the shock opening angle, 
as a function of the asymptotic velocity, computed in radians at 
r = 4.43 M and t = 10000 M, Shown with different lines are 
the values for a Schwarzschild black hole (red solid line) and for a 
Kerr black hole (blue dashed line). The points with Voo =0.3 and 
Voo = 0.4 correspond to models which are flip-flop unstable and 
they are located close to the local maximum of the sh ock opening 
angle. This results confirms what already reported by iLivio et al .1 
(1 199 lh . namely that the instability manifests when the shock open- 
ing angle is larger than a given threshold. Once the instability is 
fully developed, the mass accretion rate experiences high amplitude 
oscillations (for any value of the spin considered). This is shown 
in the bottom panel of Fig. [4] reporting the evolution of the one- 
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2000 4000 6000 8000 10 4 

t/M 

Figure 4. Top panel: shock opening angle in radians as computed at 
r = 4 ASM, t = 10000 M, as a function of Voo. Shown with different 
lines are the values for a Schwarzschild black hole (red solid line) and for 
a Kerr black hole (blue dashed line). Bottom panel: time evolution of the 
mass accretion rate, computed at r = 6.08Af, for different value of Voo 
and a rotating black hole with spin a/M = 0.9. The mass accretion rate 
manifests high amplitude oscillation in the two flip-flop unstable models 
with Voo ~ 0.3 - 0.4. 



dimensional mass accretion rate M defined as 



M 



ot^/ipW V 



11 
a 



(14) 



for different choices of Voo ; the specific values reported in the fig- 
ure refer to a shell at r = 6.08 M and a rotating black hole with 
spin a/M = 0.9. It is therefore evident that for small values of 
Voo, both the mass accretion rate and the shock opening angle are 
increasing functions of Voo ■ However, the accretion is no longer ef- 
ficient above a critical value of Voo — 0.4, and any further increase 
of Voo causes both the accretion rate and the maximum rest-mass 
density in the shock cone to decrease. 



3.2 QPOs in the shock cone 

As mentioned in the Introduction, one of the most important results 
presented in this work is about the development of QPOs in the 
shock cone that forms in the downstream region once the system 
has relaxed to the stationary state. Note that because a stationary 
state is necessary for the development of the QPOs, the latter are not 
found when the flip-flop instability develops. Although the occur- 
rence of QPOs in Bondi-Hoyle accretion has not bee n discussed be- 
fore in previous numerical investigations, Fig. 8 o f iFont & Ib anez 
dl998ij) does shows a clear periodic behavior of the accretion rate. 
So, although we can claim to be the first ones to have pointed out 
the appearance of this effect and to have investigated it in detail, we 
cannot claim to be the first to have shown found it from numerical 
simulations. 

In order to extract as much information as possible about the 
emerged phenomenology and provide a first physical explanation, 
we have carried out an extensive Fourier analysis of the differ- 
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Figure 5. Power spectra of the timeseries of the rest-mass density as computed at the center of shock cone for a flow with asymptotic velocity Voo = 0.1. 
The left and right panels refer to black holes with spins for a/M = 0.0 and a/M = 0.9, respectively, while the mass is the same in the two cases and 
M = 10 A/q. The different curves refer to timeseries recorded at different radial positions but have the same azimuthal position. 



ent dynamical quantities. The two panels of Fig. [5] for instance, 
show the power spectra obtained from the timeseries of the rest- 
mass density as computed in a region which is in the middle of the 
shock cone for a flow with asymptotic velocity Voo =0.1. The left 
and right panels refer to black holes with spins for a/M = 0.0 
and a/M — 0.9, respectively, while the mass is the same in the 
two cases and equal to M = 10 A/ Q . The different curves refer to 
timeseries recorded at different radial positions but have the same 
azimuthal position. After performing a series of tests at different 
resolutions we have estimated the error bar in the measure of each 
frequency to be ~ 3Hz. 

A number of comments should be made regarding the spectra 
reported in Fig. [5] The first one is that the modes of oscillation are 
essentially independent of the radial position, i.e., the power spectra 
at different radii overlap extremely well. This is a clear indication 
that the modes are local waves but they are global eigenmodes of 
the system. 

The second comment is that, among the modes reported in 
Fig. [5] some are genuine eigenmodes, while others are simply the 
result of nonlinear couplings. In particular, the modes with frequen- 
cies /i —21 Hz and /a = 34 Hz in the left panel of Fig.|5]are gen- 
uine eigenmodes, while those at 55, 68, 83 and 106 Hz are given, 
within a few percent error, by nonlinear couplings of the eigen- 
modes, i. e., f\ + /2, 2/2, 3/2 — fi, and 3/2 , respectively. Similarly, 
the modes /1 = 40 Hz and = 64 Hz in the right panel of Fig. [5] 
are genuine eigenmodes, while the modes at 81, 104 and 124 Hz 
are given, within few percent errors, by 2/i, fi + fa, and 3/i, 
respectively. We recall that this behavior is typical of physical sys- 
tems g overned by non linear equ ations in the limit of small oscil- 
lations dLandau & Lifshitjl 1976h and has already been pointed out 
bv lZanotti et alj J2005h in the context of oscillation modes in thick 
accretion discs around black holes. It should also be remarked that, 
while a large number of modes due to nonlinear coupling exists, the 
identification becomes difficult for frequencies larger than 100 Hz. 
Similarly, the low-frequency mode at ~ 5 Hz, which is visible in 
both the panels of Fig. [5] could be a genuine mode but it appears 
with much smaller intensity and not all radii. As a result, its classi- 
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Figure 6. Spatial distribution of the power spectral density at / = 124Hz 
for a black hole of mass M = 10 Mo and spin a/M = 0.9 (cf. right 
panel of Fig. [3). A very similar behaviour would be shown in the case of a 
nonrotating black hole. 



fication as a genuine mode will require more accurate simulations 
and on much longer timescales. 

The global nature of the oscillation modes inside the shock 
cone is clearly shown in Fig. [6] which represents the power spectral 
intensity of the most powerful mode in the right panel of Fig. [5] 
namely of the mode at / = 124 Hz (A very similar behaviour can 
be seen in the case of a nonrotating black hole.). In other words, for 
each grid-cell inside the shock cone we store the evolution of the 
rest-mass density and compute from this the power spectral density. 
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We then consider a single frequency and study how its intensity is 
distributed (in arbitrary units) inside the shock cone. 

Note that the intensity is computed over the whole compu- 
tational domain, but it has a non-negligible amplitude only inside 
the shock cone, and this amplitude becomes increasingly stronger 
near the black hole (We recall that in general we do not expect the 
amplitude to be constant in the shock cone, as this depends on the 
specific eigenfunction of the mode). What shown in Fig. [6] is not 
specific to the mode at / = 124 Hz and we have verified that the 
spatial distribution of the intensity manifests a similar pattern for 
all of the other modes. Besides providing evidence that these are 
indeed global oscillation modes of the flow, Fig. [6] also highlights 
that QPOs can be excited in the shock cone of a Bondi-Hoyle type 
accretion and that these are particularly stronger close to the black 
hole. 

Interestingly, when the flip-flop instability is triggered and de- 
velops, the power spectra inside the shock cone change consid- 
erably and in these cases only the periodicity of the shock-cone 
oscillation can be found, which is then accompanied by the usual 
nonlinear couplings. In general, therefore, the effect of the flip-flop 
instability is that of suppressing most of the internal modes of os- 
cillations. A more detailed analysis of this process, as well as of the 
flip-flop instability will be presented in a future work. 

A linear perturbative analysis of the stationary pattern would 
certainly be of great interest for three important reasons. Firstly, in 
order to identify unambiguously which of the modes are genuine 
(or fundamental) and which are the effect of non-linear coupling. 
Secondly, to compute the expected eigenfrequencies and eigen- 
functions of the oscillations we have discovered numerically. Fi- 
nally to clarify why some modes are more efficiently excited than 
others. However, such perturbative analysis is probably going to be 
extremely complicated, mainly because of the absence of an an- 
alytic solution for the stationary flow inside the shock cone. To 
counter these limitations and perform a first rudimentary analysis 
of the eigenfunctions, we report in Fig.|7]a map of what can be re- 
garded as a numerical eigenfunction for the model C2 (Voo =0.1 
and spin a/M — 0.0). In fact, the two panels of Fig.Qshow the dif- 
ference p(r, <t>)—po(r, cj>), where po(r, <j>) is the density distribution 
when the flow is stationary {i.e., for t > 1500 M). Considering, for 
example, the left panel, which corresponds to t = t\ = 1871M 
and using Cartesian coordinates, we show in the upper part of the 
figure, i.e., for y/M > 0, the density perturbation along the shock 
is negative close to the black hole, while it is positive far from it. A 
specular behavior is detected in the lower part of the figure where 
y/M < 0. The computed perturbation oscillates in time with a 
given period P, and the right panel of the figure shows a snapshot 
attimet2 = ti+P/2 = 2162M, i.e., at half a period of oscillation. 
Although heuristic, this procedure allows to compute the period of 
oscillation of the perturbation with a very good accuracy. The cor- 
responding frequency, again computed assuming M = IOMq, is, 
in fact, / ~ 35 Hz and it matches very well with the most powerful 
mode, that with frequency f%, reported in the left panel of Fig. [5] To 
highlight the dependence of the eigenfrequencies on the black hole 
spin, we report in Table[2]fheir values as measured within the shock 
cone from the evolution of the rest-mass density. The mass of the 
black hole is assumed to be M = 10 Mq, but the frequencies can 
be easily computed for an arbitrary mass since they scale linearly 
with the black hole mass. The data in Table[2]suggests that a linear 
scaling is present also for the black hole spin, but clearly more data 
is necessary to confirm this scaling, which we will investigate in a 
future work. 

As a final remark we note that, overall, the phenomenology 



Table 2. Frequencies of the first two genuine modes measured within the 
shock cone from the evolution of the rest-mass density. The values refer 
to different black hole spins and different asymptotic velocities, while the 
mass of the black hole is always assumed to be M = 10 Mq . Note that the 
frequencies scale linearly with M and the values reported have an error bar 
of ~ 3 Hz. 



a/M 


Vac 


fl [Hz] 


h [Hz] 


0.0 


0.1 


21 


34 


0.5 


0.1 


32 


53 


0.9 


0.1 


10 


64 


0.0 


0.2 


32 


46 


0.9 


0.2 


16 


81 



we have found show strong similarities with that rep orted for per- 
turbe d relativistic tori orbiting aro und a black hol e I Zanotti et al.l 
120031) . In that case it was shown bv lRezzolla et al] d2003l) through 
a perturbative investigation that the eigenfunctions and eigenfre- 
quencies were those corresponding to the p modes of the torus. 
Such modes, also known as inertial-acoustic waves, are closely re- 
lated to the propagation of sound waves in the perturbed fluid and 
have pressure gradients and centrifugal forces as the main restor- 
ing forc es. The major d i fferen ce of the physical conditions investi- 
gated by Rezzolla et al. d2003l) with the ones considered here is that 
centrifugal forces play a negligible role in Bondi-Hoyle accretion 
flows. Apart from this, however, the physical nature of the modes 
is essentially the same in the two cases. When applied to the oscil- 
lations of a thick disc around a compact object, such modes moti- 
vated the proposal of a model for explaining the detection of High 
Frequency Qua si Periodic Oscillatio ns (HFQPOs) in the spectra of 
X-ray binaries dRezzolla et ai]|2003l) . 



4 ASTROPHYSICAL APPLICATIONS 
4.1 The case of Sgr A* 

A possible application of the results discussed in the previous Sec- 
tion and, in particular, of the development of QPOs in the down- 
stream part of a Bondi-Hoyle flow, is offered by the source Sagit- 
tarius A* (or Sgr A*). Such source is located at the centre of our 
Galaxy and is widely believed to be a supermassive black hole 
with an estimated mass of ~ 4.1 ± 0.6 x 10 6 Mq dGhez et all 
Gillessen et al. IIH1). Sgr A* is one of the most intensely 
observed astronomical objects and the phenomenology of its emis- 
sion is both rich and particularly complex. Inter estingly, however, 
QPOs from its e mission have b een ob served bv lAschenbach et akl 
d2004l) (see also lAschenbachl l l2009t)). when analyzin g two dis- 
tinct ne^^nfrared_anJ Xj^ray flares " fcenzel et alj|2003l) . In partic- 
ular, lAschenb^Fetal] |200J) reported five different peaks at peri- 
ods of 100 s, 219 s, 700 s, 1150 s, and 2250 s in the power spectral 
density of such flares. Soon after. lAbramowicz et al. I (l2004l) noticed 
that (1/700) : (1/1150) : (1/2250) « 3 : 2 : 1, that is, the pe- 
riods a re in a harmo nic ratio. Indeed, while discussing the same 
source, lTor6kN2005l) confirmed the measurement of a double peak 
of QPOs in Sgr A* in the same typical 3 : 2 ratio observed in sev- 
eral low mass X-ray binaries a nd sug gested the epicyclic resonance 
model by I Abramowicz et al. (2003) as a possible explanation for 
this phenomenology. Chan et all J2009b . on the other hand, per- 
formed two dimensional magnetohydrodynamical simulations of a 
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Figure 7. Numerical eigenfunctions computed as p(r, <j>) — po(r, <j>), where p~o(r, <j>) is the density distribution when the flow is stationary (i.e., for t > 
1500 M). The left and the right panels report two snapshots at t = 1871M and t = 2162M, and are separated by half a period of oscillation; both refer for 
to a black hole with spin a/M = 0.0 and Voa =0.1. 



putative accretion disc in Sgr A* and showed that the QPOs fre- 
quencies could be due to non-axisymmetric density perturbations, 
and which may be used to set a lower limit on the orbital period at 
the innermost stable circular orbit. 

In addition to QPOs, Sgr A* is likely to have a Bondi- Hoyle 
accret ion flow. This possibility was first pointed out by IMelial 
who, after looking at the broad He I, Bra and Br7 emission 
lines, inferred that there is a strong circumnuclear wind near the 
dynamical center of the Galaxy. In this scenario, the central black 
hole is assumed to be fed by stellar winds within several a rcseconds 
from the compact object dShcherbakoy & Baganofd201(jh. In order 
to investigate this possibility further, Ruffert & Melial d 1994) per- 
formed the first three-dimensional hydrodynamical simulations of 
this system and computed the line-integrated flux assuming that the 
emission is dominated by bremsstrahlung. B y collecting informa- 
tion coming from more r e cent o bservations, ICuadra et alj d2006t) 
and ICuadra & Navakshinl d2006t) performed SPH simulations of 
wind accretion onto Sgr A* , including optically thin radiative cool- 
ing and finding, among other results, that most of the accreted gas 
is hot and has a nearly-stationary accretion rate. Although the ob- 
served luminosity is some orders of magnitude smaller than what 
is predicted by the Bondi-Hoyle model, this may be due to a low 
radiative efficiency of the accretion flow, rather than to a failure of 
the Bondi-Hoyle model itself. 

Should the two observational features mentioned above be 
confirmed by additional and more accurate observations, there 
would be a unique astronomical realization of a source for which 
both QPOs and Bondi-Hoyle accretion are simultaneously present. 
It should be stressed that the typical lengthscale at which the QPOs 
develop in our simulations is ~ 20 M, which thus corresponds to 
~ 4 x 10 -6 pc for the estimated mass of Sgr A*. This lengthscale 
is much smaller than the closest approach to the Galactic Center of 
the closest s tar SI 6, which is around ~ 2 x 10 -4 pc, as reported by 
iGhez et al] d2005h . Therefore, the Bondi-Hoyle accretion flow that 
could potentially be present in Sgr A*, would take take place in a 
gas-dominated region very close to the central black hole and far 
from the orbits of the S-stars. 

Interestingly, after analyzing recent data about a high-level X- 



ray activity of Sgr A* ob s erved with XMM-N ewton dPorquet et alj 
20081; lAschenbachll2009h . | Aschenbachl d2009l) reports nine promi- 
nent peaks, name d v\ . . . vg. Out of th em, only three are identified 
as fundamental by |Aschenbacr]d2009l) , namely ue, vi and vg, while 
the other six are obtained as linear combinations of th e fundamen- 
tal one s, like, for instance, uq + vi, uq + vg, v~i — vg. I Aschenbachl 
d2009h interprets the three fundamental modes in terms of orbital, 
radial epicyclic and vertical epicyclic frequencies. In doing so, a 
specific radius where the oscillation is excited must be provided, 
the above osscillation frequencies being intrinsecally local. On the 
contrary, no radial specification is required within our interpreta- 
tion, since the modes are global and just c onfined within the shock 
cone. Finally, althoug h i Aschenbachl d2009h does not give any expla- 
nation for the other six modes, they are naturally provided within 
our interpretation of a Bondi-Hoyle accretion in which the down- 
wind shock cone acts as a cavity, generating a whole series of non- 
linear couplings among few genuine and trapped pressure modes. 



4.2 QPOs in HMXBs 

A second application of our result is in principle represented 
by QPOs observed in the spectra of high mass X-ray binaries 
(HMXBs), which are compose d of a compact ob ject and of an early 
type (OB) star. The catalog bv lLiu et alj d2006h lists 114 HMXBs 
candidates in the Galaxy, some of whic h are believed t o contain a 
bl ack hole, like Cyg X-l, NGC 520 4 dLiu et ai]|2004h M 33 X- 
7 dPietsch et al.l2004) . M 101 ULX-1 dMukai et alj2005h . Because 
in these systems the accretion flow occurs preferentially in the form 
of an accretion wind rather than in the form of an accretion disc, the 
Bondi-Hoyle accretion flow has bee n traditionally considered very 
relevant for them (see the review bv lEdgarl 12004) for the applica- 
tion of the Bondi-Hoyle solution to binary systems). However, the 
QPOs that have been detected in HMXBs have frequencies that are 
seen to lie in the range of 1 mHz to 400 mHz, with the remarkable 
exceptio n of XTE J011 1.2 -73 17, which shows a QPO feature at 
1.27 Hz dKaur et alj|2007l) . On the other hand, the QPOs we have 
computed assuming a black hole of mass M = 10 Mq have fre- 
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quencies which lie in the range from 1 Hz to 500 Hz , and therefore 
that overlap only marginally with the observed ones. 



5 CONCLUSIONS 

Although viscous accretion discs represent the most natural chan- 
nel by means of which compact objects accrete large amounts of 
matter with significant angular momentum, non-spherical accretion 
flows appear in all those situations in which the accreting matter 
has only a modest amount of angular momentum, such as in winds. 
Bondi-Hoyle accretion is the most representative example of this 
type of flow and its study via numerical investigation has a long his- 
tory, both in Newtonian and in general-relativistic physics. We have 
reconsidered this old problem by performing new two-dimensional 
and general-relativistic simulations onto a rotating black hole. Be- 
sides recovering many of the features of this flow which were dis- 
cussed by other authors over the years, we have also pointed out a 
novel feature. More specifically, we have shown that under rather 
generic conditions a shock cone develops in the downstream region 
of the flow and that such a cone acts like a cavity trapping pressure 
modes and giving rise to QPOs. 

These modes are global in the sense that they represent har- 
monic oscillations across the cavity, but have amplitudes that are 
larger in the region very close to the black hole, i.e., for r < 10M. 
While the black hole spin influences the absolute frequencies of the 
trapped modes, which lie in the range from 1 Hz to 500 Hz for a 
representative black hole of mass M = 10 Mq, it does not affect 
the fact that they appear in a series of integer numbers 1:2:3. 
This is because nonlinear coupling of modes is common in systems 
governed by nonlinear equations, and once a mode is excited, all of 
its integer multiples are also excited, thus producing a wide range 
of possible ratios of integer numbers. 

In addition to pointing out this novel feature of the 
Bondi-Hoyle accretion, we have discussed its possible appli- 
cation to_Jte_pJienomenology reported in SgrA*, wher e both 
QPO s dAschenbach et alj[2004l) and Bondi-Hoyle accretion l lMelial 
Il992h are likely to be present. Remarkably, a large fa mily of linear 
combi nations of modes has been identified recently bv lAschenbachl 
d2009h after analyzing data about a hi gh-level X-ray activ ity of 
Sgr A* observed with XMM-Newton jporquet et alTl2008h - This 



phenomenology could indeed be interpreted in terms of a shock 
cone that behaves like a cavity and that generates a whole class of 
nonlinear coupling among pressure modes. 

Finally, we have provided the first evidence for the occurrence 
in a general relativistic context of the so called flip-flop instability 
of a Bondi-Hoyle flow. More specifically, we have shown that for 
a fixed choice of the black hole spin and of the sound speed, there 
exists a critical value of the asymptotic flow velocity at which the 
shock cone undergoes large-scale and coherent oscillations. When 
this happens, the opening angle of the shock cone reaches its max- 
imum value, the accretion rate increases considerably and is no 
longer stationary. A more comprehensive analysis is needed to clar- 
ify the physical nature of the instability and its dynamics in the 
relativistic regime. This will be the focus of a future work. 

Finally, among the future improvements of our work we plan 
to investigate the effects of a magnetic field, which is certainly 
likely to play a role but whose effective contribution has not been 
considered yet in numerical simulations of Bondi-Hoyle accretion 
flows. In particular, the potential interplay of the magnetorotational 
instability (MRI) with the flip flop instability has to taken into ac- 
count, though it is not clear whether a local instability such as the 



MRI can substantially affect a global instability. As far as the oc- 
currence of QPOs, on the other hand, we believe that it will remain 
essentially unchanged. The reason for this is that as long as a shock 
cone forms in the downstream part of the flow, QPOs will be natu- 
rally excited and trapped, the only difference being that they will be 
of magnetosonic nature rather than of a sonic one, and hence with 
eigenfrequencies that will depend not only on the properties of the 
shock-cone (and thus of the fluid) but also on the strength of the 
magnetic field. 
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